In [1]:
%matplotlib inline
import cclib
import numpy as np
from scipy import linalg
from matplotlib import pyplot as plt
import seaborn as sns
sns.set_style("white")
In [2]:
def modiagram(mol, frag1, frag2):
"""Return a molecular orbital diagram by projecting the molecular orbitals of a molecule onto the ones of two
different fragments.
Fragments should have the same basis, geometry and atomic ordering as the molecule.
Parameters
----------
mol : str
File name for the output calculation of the molecule to be read by cclib.parser.ccopen.
frag1 : str
File name for the output calculation of the first fragment to be read by cclib.parser.ccopen.
frag2 : str
File name for the output calculation of the second fragment to be read by cclib.parser.ccopen.
Returns
-------
mole : array_like
Energies the molecular orbitals of the molecule.
frag1e : array_like
Energies the molecular orbitals of the first fragment.
frag2e : array_like
Energies the molecular orbitals of the second fragment.
U : array_like
Transformation matrix from fragment molecular orbitals to molecular orbitals of the final molecule."""
pmol = cclib.parser.ccopen("H2.out")
pfrag1 = cclib.parser.ccopen("H2_f1.out")
pfrag2 = cclib.parser.ccopen("H2_f2.out")
dmol = pmol.parse()
dfrag1 = pfrag1.parse()
dfrag2 = pfrag2.parse()
# The variables below hold MO coefficients of the supermolecule and its fragments in terms of basis set
# functions.
Cm_a = dmol.mocoeffs[0]
if len(dmol.mocoeffs) > 1:
Cm_b = dmol.mocoeffs[1]
Ck_a = linalg.block_diag(dfrag1.mocoeffs[0], dfrag2.mocoeffs[0])
if len(dfrag1.mocoeffs[1]) > 1 and len(dfrag2.mocoeffs[1]) > 1:
Ck_b = linalg.block_diag(dfrag1.mocoeffs[1], dfrag2.mocoeffs[1])
# Sm below holds the basis set overlaps, taken from the molecule calculation. The overlaps taken from the
# fragment calculations would be
#
# Sk = linalg.block_diag(dfrag1.aooverlaps, dfrag2.aooverlaps)
#
# but this matrix would lack many cross-terms.
Sm = dmol.aooverlaps
# The following would produce a block matrix [[0, A], [A, 0]]. This shows we have the right basis and the
# right ordering of it.
#
# plt.matshow(np.abs(Sk - Sm))
# plt.colorbar()
# Next, a matrix will be calculated that holds the transformation from the MOs of fragments to the MOs of the
# molecule.
Cf_a = np.linalg.multi_dot([Cm_a, Sm, Ck_a.transpose()])
# if Ck_b and Cm_b are defined:
# Cf_b = np.linalg.multi_dot([Cm_b, Sm, Ck_b.transpose()])
# The following code depicts the matrix and shows that the resulting MO orbitals are normalized:
#
# plt.matshow(Cf_a)
# plt.colorbar()
# print(np.sum(Cf_a**2, axis=0))
return dmol.moenergies, dfrag1.moenergies, dfrag2.moenergies, Cf_a
In [3]:
def plot_modiagram(mole, frag1e, frag2e, U, thresh=.0817, mowidth=1., mosep=.5, cog=0.):
"""Plot a molecular orbital diagram.
Parameters
----------
thresh : float
Minimum value for drawing dashed lines between fragment and molecule molecular orbitals.
mowidth : float
Width in the x-axis of each plotted molecular orbital.
mosep : float
Distance in the x-axis of minimum separation between molecular orbitals.
cog : float
x position for the center of the graph."""
# We check first multiplicities.
# https://docs.scipy.org/doc/numpy-1.10.1/reference/generated/numpy.unique.html
# umos = [mos[0]]
# dmos = default{0: 1}
# last_e = mos[0]
# for e in mos[1:]:
# if e == last_e:
# if
def _plot_levels(levels, mowidth=mowidth, mosep=mosep, cog=cog):
for e in levels:
# Write a line at energy e with width mowidth
plt.plot((cog - (mowidth - mosep)/2., cog + (mowidth - mosep)/2.), (e, e), 'k-')
_plot_levels(mole, cog=cog)
_plot_levels(frag1e, cog=cog - mowidth)
_plot_levels(frag2e, cog=cog + mowidth)
frag1n = len(frag1e)
frag2n = len(frag2e)
for i, e in enumerate(mole):
for j in range(len(U[i, :frag1n])):
csqr = U[i, j]**2
if csqr > thresh:
begin = np.array([cog - mowidth + (mowidth - mosep)/2., frag1e[j]])
end = np.array([cog - (mowidth - mosep)/2., e])
plt.plot((begin[0], end[0]), (begin[1], end[1]), 'k--')
plt.annotate('{:.0f}%'.format(csqr*100.), xy=(begin + end)/2.)
for j in range(len(U[i, frag1n:(frag1n + frag2n)])):
csqr = U[i, j]**2
if csqr > thresh:
begin = np.array([cog + (mowidth - mosep)/2., e])
end = np.array([cog + mowidth - (mowidth - mosep)/2., frag1e[j]])
plt.plot((begin[0], end[0]), (begin[1], end[1]), 'k--')
plt.annotate('{:.0f}%'.format(csqr*100.), xy=(begin + end)/2.)
In [4]:
ethresh = 5.
mole, frag1e, frag2e, U = modiagram("H2.out", "H2_f1.out", "H2_f2.out")
mole_a = mole[0][mole[0] < ethresh]
frag1e_a = frag1e[0][frag1e[0] < ethresh]
frag2e_a = frag2e[0][frag2e[0] < ethresh]
plot_modiagram(mole_a, frag1e_a, frag2e_a, U)
plt.ylabel("Energy (eV)")
plt.xticks([])
Out[4]: